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NOTATION 


A  =  Vehicle  cross  sectional  area 

a  s  Semi-major  axis 

B*  =  Eclipse  entrance/exit  angle 

Qi  s  Drag  coefficient 

D  s  Drag  force 

E  =  Specific  mechanical  energy 

F  s  Applied  force  on  vehicle 

fe  =  Fraction  of  orbit  in  eclipse 

fs  =  Fraction  of  orbit  in  sun 

gf.  s  Pound  mass/pound  weight  conversion  factor 

go  s  Gravitational  acceleration  at  the  earth’s  surface 

h  =  Vehicle  altitude 

i  =  Inclination 

Isp  5  Specific  impulse 

Ai  2  Inclination  change  imparted  to  the  orbit 

h  s  C20  gravitational  coefficient  (oblateness) 

m  2  Vehicle  mass 

m  =  Propellant  mass  flow  rate 

Mf  =  Dry  vehicle  mass 

M[  2  Wet  vehicle  mass 

Mt  2  Thermal  management  system  mass 

P  2  Available  power 

Pi  2  Initial  system  power 

Pr  2  Period  of  transfer  orbit 

q  2  Dynamic  pressure 

r  2  Instantaneous  radius  of  spacecraft  orbit 

r  =  Rate  of  change  of  radius 

ri  2  Radius  of  initial  orbit 
X2  =  Radius  of  final  orbit 
Re  2  Radius  of  Earth 
T  2  Total  thrust 

t  =  Time  since  beginning  of  transfer 

Teff  =  Effective  thrust 


Tn  =  Normal  thrust 

to  2  Time  of  vernal  equinox  (March  21) 

Tp  2  In-plane  thrust 

V  2  Orbit  velocity 

vt  2  Initial  orbit  velocity 
V2  =  Final  orbit  velocity 
AV  2  Velocity  change  imparted  to  the  orbit 
Wt  2  Thermal  managment  system  specific  mass 
a  2  Out-of-plane  steering  angle 
as  2  Solar  right  ascension 
6  2  Solar  declination 

e  2  Obliquity  of  ecliptic 

Y  s  Flight  path  angle 

T[ppU  2  Power  processor  efficiency 

T],  2  Thruster  efficiency 

\i  2  Earth’s  gravitational  parameter 

p  2  Atmospheric  density 
• 

0  2  Angular  rate  of  vehicle 

CDs  2  Angular  motion  of  sun 

O  2  Right  ascension  of  ascending  node 

• 

=  Rate  of  change  of  £2 
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INTRODUCTION 


The  increasing  technological  maturity  of  high  power  (>20  kW)  electric  propulsion  devices  has  led  to 
renewed  interest  in  their  use  as  a  means  of  efficiently  transferring  payloads  between  earth  orbits.  Several 
systems  and  architecture  studies  have  identified  the  potential  cost  benefits  of  high  performance  Electric 
Orbital  Transfer  Vehicles  (EOTVs)1*2.  These  studies  led  to  the  initiation  of  the  Electric  Insertion  Transfer 
Experiment  (ELITE)  in  1988.3  Managed  by  the  Astronautics  Laboratory,  ELITE  is  a  flight  experiment 
designed  to  sufficiently  demonstrate  key  technologies  and  options  to  pave  the  way  for  the  full-scale 
development  of  an  operational  EOTV. 

An  important  consideration  in  the  development  of  the  ELITE  program  is  the  capability  of  available 
analytical  tools  to  simulate  the  orbital  mechanics  of  a  low  thrust,  electric  propulsion  transfer  vehicle. 
These  tools  are  necessary  not  only  for  ELITE  mission  planning  exercises  but  also  for  continued,  efficient, 
accurate  evaluation  of  DoD  space  transportation  architectures  which  include  EOTVs.  This  paper  presents 
such  a  tool:  the  Electric  Vehicle  Analyzer  (EVA). 

BACKGROUND 

To  properly  understand  the  need  for  EVA  and  its  modeling  capabilities,  it  is  necessary  to  look  at  some 
background  information  which  has  shaped  its  development. 

EOTV  Design  and  Operation 

An  Orbital  Transfer  Vehicle  (OTV)  is  a  vehicle  used  to  transfer  payload  (e.g.,  a  satellite)  from  a  low 
earth  orbit  into  a  higher  orbit.  One  of  the  most  common  destination  orbits,  as  well  as  one  of  the  most 
demanding  on  the  OTV,  is  the  geosynchronous  orbit.  A  circular  orbit  in  the  plane  of  the  equator  at  an 
altitude  of  35000  km,  with  a  period  of  exactly  one  day.  This  orbit  is  popular  for  communications 
satellites. 

The  main  difference  between  conventional  OTVs  and  EOTVs  is  the  type  of  propulsion.  Conventional 
OTVs  use  chemical  propulsion,  whereas  EOTVs  utilize  externally  supplied  electrical  power  instead  of  a 
chemical  reaction.  There  are  different  options  possible  to  supply  the  electric  power  including  photovoltaic, 
solar  dynamic,  and  nuclear  power  systems.  ELITE  has  baselined  photovoltaic  solar  arrays  to  provide 
electrical  power. 

pm, 

W 

Figure  1 

OTV  Comparison 
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Although  the  theoretical  Isp  limits  on  chemical  propulsion  (=500  sec)  do  not  apply  to  electric 
propulsion,  allowing  IspS  up  to  5,000  seconds,  photovoltaic  power  limitations  restrict  thrust  levels  to  less 
than  10N  (=2  lb).  Short,  impulsive  bums  which  are  characteristic  of  chemical  propulsion  are  not  possible 
at  such  low  thrust  levels,  and  continuous  thrusting  throughout  the  mission  becomes  the  standard  mode  of 
operation.  Hohmann  transfers,  which  consist  of  two  impulsive,  high-thrust  bums  (Figure  2),  give  way  to 
trajectories  that  slowly  spiral  upward  to  the  destination  orbit,  circling  the  earth  hundreds  or  even  thousands 
of  times  along  the  way. 


Figure  2 

Hohmann  vs  Spiral  Trajectory 

Desired  Analytical  Capability 

In  order  to  do  their  job,  EOTV  mission  planners  will  have  to  have  knowledge  of  the  effects  of 
operating  a  vehicle  in  such  a  slow,  low-thrust,  spiral  trajectory.  Some  of  these  effects  include  gravity 
losses,  atmospheric  drag,  solar  occultation,  and  radiation  damage  to  the  solar  arrays. 

Gravity  Loss.  One  consequence  of  flying  a  low-thrust  spiral  trajectory  rather  than  the  conventional 
two  impulse  bum  Hohmann  transfer  is  gravity  loss.  In  the  ideal  Hohmann  scenario,  an  impulse  bum 
takes  place  at  periapsis  and  at  apoapsis  where  the  thrust  vector  is  perpendicular  to  the  gravity  vector  of  the 
central  body.  All  thrust  applied  can  be  used  to  add  energy  to  the  orbit.  Longer  bum  times  create  the 
situation  shown  in  Figure  3;  in  which  the  bum  takes  place  through  points  on  the  orbit,  at  which  the  flight 
path  angle  is  not  zero.  A  portion  of  the  thrust  must  therefore  by  used  to  overcome  the  local  effect  of 
gravity.  This  manifests  itself  as  a  loss,  and  for  low  earth  to  geosynchronous  missions  can  increase  total 
velocity  (Av )  required  to  complete  the  mission  by  30%  . 
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Figure  3 

Impulsive  vs  Non-Impulsive  Burn 


Atmospheric  Drag.  With  current  state  of  the  art  solar  array  technology,  array  areas  on  the  order  of 
hundreds  of  square  meters  are  necessary  to  provide  enough  electric  power  to  produce  reasonable  thrust 
and  trip  times  for  EOTVs.  These  large  areas,  combined  with  the  low  thrust  of  the  electric  thrusters,  make 
atmospheric  drag  at  low  earth  orbits  a  significant  factor  in  mission  planning.  These  factors  may  even 
dictate  minimum  deployment  altitudes  below  which  the  thrust  to  drag  ratio  will  be  unacceptable,  causing 
decay  of  the  vehicle’s  orbit  and  ultimately  atmospheric  reentry. 

Occultation.  Time  spent  in  the  Earth’s  shadow,  or  solar  occultation,  is  another  significant  factor  to 
consider  in  the  mission  planning  of  solar-powered  EOTVs.  Since  such  an  EOTV  derives  its  electrical 
power  from  the  sun,  it  must  coast  through  the  occultation  periods.  This  not  only  alters  the  shape  of  the 
spiral  trajectory  but  also  affects  the  trip  time.  To  illustrate  this  effect  consider  an  EOTV  starting  out  in  a 
circular  orbit  (Figure  4).  It  is  thrusting  when  in  sun;  not  thrusting  when  in  shadow.  This  period  of 


Figure  4 

Effect  of  Occultation  on  Trajectory 
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coasting  will  cause  the  orbit  to  become  slightly  elliptical.  Also,  since  the  shadow  cylinder  rotates  in  inertial 
space  as  the  Earth  moves  in  its  orbit  about  the  sun,  the  perigee  location  of  the  now  slightly  elliptical  orbit 
will  also  rotate  in  inertial  space.  Thus  as  the  EOTV  spirals  upward  toward  its  destination  orbit,  its 
trajectory  will  be  a  slightly  elliptical  spiral  rather  than  a  perfectly  circular  spiral. 

A  typical  eclipse  history  is  shown  in  Figure  5  for  a  transfer  between  low  earth  orbit  and 
geosynchronous  orbit.  Over  the  course  of  the  entire  transfer,  the  vehicle  will  experience  several  hundred 
occultation  cycles  resulting  in  a  total  trip  time  penalty  of  up  to  several  weeks. 


Cumulative 

Ecllpsa 

Duration 

(days) 


Launch  LEO 
June  19  28.5° 


Altitude  (kmxlOOO) 


GEO  Arrive 
0'  Sept  17 


Figure  5 

Typical  Eclipse  History  for  a  LEO  to  GEO  Transfer 


One  option  available  to  overcome  this  phenomenon  is  to  carry  enough  battery  power  to  continue 
thrusting  through  the  occultation  periods.  Conventional  batteries  could  provide  the  power,  however,  their 
weight  would  be  prohibitive. 


Another  factor  to  consider  is  the  combination  of  atmospheric  drag  and  occultation.  The  vehicle  will 
undergo  occultation  even  at  low  earth  orbit  where  atmospheric  drag  is  most  severe.  Depending  on  the  size 
of  the  arrays,  the  amount  of  thrust  generated  by  the  thrusters,  and  the  actual  altitude  of  the  low  earth  orbit, 
the  vehicle’s  orbit  could  easily  deteriorate  to  reentry  because  of  drag  during  an  occultation  period. 

Radiation.  At  mid  altitudes  (1,000  -  10,000  km)  lie  the  Van  Allen  radiation  belts.  Consisting  of 
trapped,  high-energy  electrons  and  protons,  these  regions  can  cause  serious  damage  to  photovoltaic  solar 
arrays,  and  vehicle  electronics.  Because  of  the  low  thrust  nature  of  EOTVs,  the  vehicle  will  spend  many 
days  in  the  Van  Allen  belts.  The  primary  effect  on  the  EOTV  itself  is  the  damage  to  solar  arrays  which 
manifests  itself  in  the  form  of  reduced  power  producing  ability.  As  the  power  degrades,  the  thrusters  will 
have  to  be  throttled  or  shut  down,  which,  in  turn,  will  have  an  effect  on  the  shape  of  the  transfer  trajectory 
and  total  transfer  time.  Shielding  can  be  added  to  the  solar  arrays  to  minimize  the  radiation  damage; 
however,  a  weight  penalty  is  incurred. 
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Each  of  these  effects,  while  not  necessarily  important  for  conventional  chemical  OTVs,  are  very 
important  to  consider  in  the  mission  planning  of  electric  OTVs.  Any  useful  mission  planning  tool  should 
address  these  factors  to  some  degree. 

Also  of  importance,  especially  for  the  ELITE  mission  planning  tool,  is  the  capability  to  model 
different  electric  propulsion  systems  and  mission  parameters.  Preferably  this  would  be  accomplished  by 
making  quick  and  easy  adjustments  to  an  input  data  file.  While  very  accurate  evaluation  of  the  various 
low-thrust  effects  presented  above  is  desirable,  a  rigorous  integration  of  complex  orbit  and  attitude 
equations  would  slow  the  program  down  so  much  as  to  be  undesirable  for  use  in  the  PC  environment. 
The  goal  was  to  have  a  code  that  would  be  as  accurate  as  possible  while  maintaining  run  times  of  less  than 
30  seconds  in  a  VAX  environment  and  also  be  PC  compatible. 

Existing  Tools 

A  survey  of  low  thrust  transfer  codes  currently  in  use  within  government  agencies  was  conducted 
and  two  validated  and  well  proven  programs  were  identified:  SPACEDRIVE  and  the  Solar  Electric 
Control  Knob  Setting  Program  for  Optimal  Trajectories  (SECKSPOT). 

SPACEDRIVE.  An  MS  DOS-based  program  developed  by  Electric  Propulsion  Lab,  Inc.  of 
Lancaster,  CA;  is  continuing  to  evolve  as  a  versatile  mission  analysis  and  propulsion  subsystem  design 
tool.4  In  an  effort  funded  by  the  Strategic  Defense  Initiative  and  managed  by  NASA/Lewis  Research 
Center,  the  Electric  Propulsion  Lab  developed  this  user-interactive  software.  The  mission  analysis  routine 
in  SPACEDRIVE  was  investigated  for  applicability  to  specific  AL  needs.  Some  simplifications  employed 
by  SPACEDRIVE  such  as  the  restriction  to  coplanar  transfers,  the  absence  of  atmospheric  drag,  and  the 
inability  to  model  solar  power  degradation  were  seen  as  detractions  to  the  tool’s  strong  point  of  user 
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friendliness.  The  particular  iterative  orbit  transfer  scheme  used  by  SPACEDRIVE  also  made  run  times 
longer  than  was  desired.  For  these  reasons,  SPACEDRIVE  was  not  used  for  our  orbital  transfer  mission 
planning  program.  However,  we  continue  to  use  SPACEDRIVE  for  its  extensive  database  and  subsystem 
design  model. 

SECKSPOT.  This  program  represents  the  other  side  of  the  analytical  spectrum  to  SPACEDRIVE. 
Written  in  the  early  1970’s  by  MIT’s  Charles  Stark  Draper  Laboratory,  SECKSPOT  optimizes  steering 
control  algorithms  to  minimize  time  of  flight  for  low  thrust  propulsion  systems.5  SECKSPOT  continu¬ 
ously  updates  the  vehicle’s  orbital  state  and  attitude  in  its  analysis  of  the  trajectory.  Among  other  things, 
the  program  takes  into  account  occultation,  power  system  degradation,  and  earth  oblateness.  It  represents 
one  of  the  highest  fidelity  low  thrust  transfer  programs  available.  The  large  number  of  differential  and 
integral  equations  which  must  be  solved,  however,  make  it  extremely  unwieldy  in  a  PC  or  even  a  VAX 
environment.  SECKSPOT,  for  this  reason,  was  determined  to  be  inappropriate  for  the  desired  capability. 

Analytical  Gap.  SPACEDRIVE  and  SECKSPOT  represent  the  only  two  well  documented  and 
proven  software  tools  available  within  government  agencies.  The  mission  planning  capabilities  of  these 
programs  represent  two  ends  of  the  analytical  spectrum.  Unfortunately,  the  needs  of  the  Astronautics 
Laboratory  fall  in  the  center  of  this  “analytical  gap”.  As  such,  it  was  determined  that  in-house  propulsion, 
orbital  mechanics,  and  software  development  expertise  should  be  pooled  to  develop  a  tool  for  parametric 
mission  studies  which  encompassed  all  of  the  capabilities  discussed  above.  This  program  came  to  be 
known  as  the  Electric  Vehicle  Analyzer. 

ELECTRIC  VEHICLE  ANALYZER  PROGRAM 


Program  Architecture 

EVA  was  designed  as  a  user-friendly,  preliminary  mission  planning  tool  that  would  be  able  to  model 
the  effects  of  a  low-thrust  trajectory  as  accurately  as  possible  while  keeping  run  times  on  the  order  of 
seconds  on  a  VAX  and  minutes  on  a  PC.  It  was  also  designed  to  have  the  ability  to  model  a  variety  of 
missions  and  types  of  electric  propulsion  systems.  It  was  structured  so  that  it  could  be  easily  upgraded  to 
include  more  user-friendly  features  or  more  accurate  mathematical  models. 

To  accomplish  the  user-friendliness  goal,  the  input  data  for  EVA  was  separated  and  placed  into  a  self- 
contained  input  file  in  namelist  format  Any  or  all  of  the  parameters  could  be  easily  changed  in  this  file  and 
then  the  program  rerun.  This  method  was  chosen  so  as  to  maximize  the  number  of  input  parameters  that 
could  be  changed  while  at  the  same  time  eliminating  the  need  to  re-input  the  entire  data  set  in  a  lengthy 
series  of  interactive  questions  at  the  beginning  of  each  run.  All  of  the  mission  parameters,  electric 
propulsion  system  specifications,  and  program  flags  are  contained  in  this  input  file.  See  Appendix  A  for  a 
complete  list  of  variables  with  definitions. 

The  main  code  is  separated  into  several  separate  subroutines.  The  workhorse  of  the  program  is 
Subroutine  Transfer  which  handles  all  of  the  orbital  transfer  calculations  (Appendix  A).  The  atmosphere 
model,  the  radiation  tables,  and  the  power  degradation  calculations  are  further  separated  into  individual 
subroutines.  The  main  reason  for  choosing  this  construction  is  to  facilitate  future  upgrades  in  the 
individual  mathematical  models  so  they  can  be  easily  incorporated  into  the  code. 
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Mathematical  Mode! 

The  primary  goal  of  the  EVA  orbital  mechanics  mathematical  model  was  to  describe  the  effects  of  low 
thrust  trajectories  as  accurately  as  possible  without  taking  a  lot  of  CPU  time.  In  order  to  accomplish  this, 
some  assumptions  were  made  about  the  nature  of  the  transfer  process.  These  will  be  discussed  along  with 
the  actual  mathematical  derivations  in  the  following  sections. 


Altitude  Raising.  At  the  heart  of  EVA  is  the  equation  which  determines  the  rate  of  change  of  the 
orbital  radius.  The  effects  of  plane  changing,  atmospheric  drag,  power  degradation,  and  other  important 
perturbations  hinge  on  this  key  equation.  The  assumption  made  here  is  that  the  vehicle  is  on  a  perfectly 
circular  spiral  trajectory.  That  is  to  say  that  at  every  altitude  along  the  spiral,  the  vehicle  trajectory  is 
tangent  to  an  equivalent  circular  orbit  of  that  altitude.  As  mentioned  earlier,  this  may  not  be  completely 
accurate  due  to  the  occupation  phenomenon,  however,  it  has  been  estimated  that  the  cumulative  effect 

during  a  low  earth  to  geosynchronous  transfer  would  result  in  an  elliptic 
geosynchronous  orbit  with  an  eccentricity  of  less  than  0.2.  Following  the  ap¬ 
proach  by  Shepard6,  the  radial  rate  equation  can  be  derived  as  follows.  The 
specific  mechanical  energy  of  the  vehicle  is  defined  as 


where  ji  is  the  earth’s  gravitational  parameter  and  a  is  the  orbit  semi-major  axis. 
Based  on  the  assumption  that  the  spiral  transfer  orbit  remains  nearly  circular  at  any 
given  instant,  the  semi-major  axis  can  be  replaced  by  circular  orbit  radius,  r,  and  the  equation  differentiated 
to  give 


dE 

dt 


M_da 
2  dt 


2a 


f^dr 
2  dt 


2r 


(2) 


For  tangentially  applied  forces,  the  rate  of  change  of  a  vehicle’s  orbital  energy  is 

dE  T7W 
m-ar=g<FV 


(3) 


where  m  is  the  mass  of  the  vehicle,  F  is  the  tangential  force,  V  is  the  orbit  velocity  and  gc  is  a  conversion 


factor.  Combining  Equations  2  and  3  gives, 

,  l  2\ 

dr 


s-^l- 


m 


For  a  near  circular  orbit,  the  velocity  can  be  approximated  by 


(4) 


(5) 


Substituting  Equation  5  into  Equation  4  and  rearranging  gives  the  instantaneous  rate  of  change  of 
transfer  orbit  radius  as  a  function  of  altitude  and  tangential  acceleration.  This  is  the  key  equation  used  for 
computing  the  spiral  orbit  altitude  increase  with  time. 


(6) 
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Another  useful  piece  of  information  regarding  the  transfer  orbit  is  the  instantaneous  flight  path  angle. 
This  can  be  directly  derived  from  Equation  4. 

2'  h 

(7) 


Y= 


*/ 

/A 


r  F 


For  constant  accelerations,  Equation  6  could  be  integrated  directly  to  determine  orbit  radius  as  a 
function  of  time.  Unfortunately,  Equation  6  becomes  non-linear  when  atmospheric  drag,  power  system 
degradation,  and  thrust  discontinuity  due  to  eclipse  are  considered.  Piecewise  integration  is  used  to  solve 
Equation  6  over  intervals  in  which  the  forces  can  be  considered  constant.  The  average  altitude  raising 
force  over  the  course  of  the  orbit  is 


F  =  TPfs-D  (8) 

where/5  is  the  fraction  of  the  orbit  spent  in  sunlight  with  thrusters  operating  and  Tp  is  the  in-plane  thrust 
component.  F  is  comprised  of  thrust  forces  and  (at  low  altitudes)  drag  forces — other  perturbation  forces 
are  neglected.  The  calculation  of  /$  is  discussed  in  the  section  on  power  system  interactions.  The 
computed  value  for  dr/dt  is  then  multiplied  by  the  orbital  period  to  find  the  change  in  altitude. 


Orbital  Period  Adjustment.  At  low  earth  altitudes  where  the  time  rate  of  change  of  the  spiral 
radius  is  small,  the  period  is  very  close  to  the  period  of  a  circular  orbit  at  that  altitude.  For  larger  values  of 
dr/dt,  however,  the  period  will  be  somewhat  longer,  and  a  correction  factor  is  introduced  so  that  the 
transfer  period  can  be  more  accurately  determined.  The  derivation  of  the  orbital  period  equation  follows 
from  Figure  7.  Approximating  to  near  circular  orbits  yields 


(9) 


The  instantaneous  radius  during  an  orbit  revolution  can  be  approximated 
r  =  rx  +  rt 


(10) 


where  r  is  the  average  calculated  rate  of  change  given  by  Equation  6  and  ri  is  the  radius  at  the  start  of  the 
orbit.  Combining  Equations  9  and  10  yields 


r,+  r  t 


(11) 
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Figure  7 


The  Spiral  Orbit  Trajectory 
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Inclination  Changes.  Since  most  satellites  start  out  in  a  low  earth  orbit  that  is  inclined  to  some 
degree,  and  since  many  orbital  transfer  missions  have  a  destination  of 
geosynchronous  orbit  at  0°  inclination,  any  mission  planning  tool  should  have  the 
ability  to  model  inclination  changes.  This  requirement  introduces  the  need  to  use 
some  thrust  in  an  out-of-plane  direction.  This  is  readily  handled  by  defining  the 
desired  thrust  vector  out-of-plane  steering  angle  a.  The  total  thrust  available  from 
the  electric  thruster  devices  can  be  computed  knowing  the  available  power,  system 
efficiencies,  and  specific  impulse. 


m  = 


2TltTl 


PEL 


(go^sp) 


(17) 

(18) 


T  -  Isp  go  m 

The  normal  thrust  used  for  inclination  changing,  therefore,  is  given  by 

Tn=T  sina  (19) 

where  T  is  the  total  delivered  thrust.  While  the  available  in-plane  (altitude  raising)  thrust  is  simply 
TP  =  Tcosa  (20) 

Using  the  relationship  developed  by  Edelbaum7  for  the  low  thrust  delta-velocity  requirement  between  two 
circular,  non-coplanar  orbits, _ 

AV  =  /\J  V12+V22-2V1V2cos^Ai)  (21) 

the  incremental  inclination  change  during  an  orbit  revolution  can  be  determined  by  solving  for  Ai. 
2  -u  v » 1  T  »2 


Ai  =  —  cos 
n 


i(W>V22-AV2 


2V,V2 


(22) 


The  velocities  at  the  start  and  end  of  the  orbit  revolution  are  determined  from  Equation  5  as  discussed  in 
the  previous  section.  The  aV  can  be  calculated  from  the  rocket  equation 

AV  =  'sPS.^)  (23) 

where  the  difference  between  Mi  and  Mp  is  the  amount  of  propellant  consumed  during  the  last  orbit  which 
can  be  calculated  knowing  the  characteristics  of  the  propulsion  system.  It  is  important  to  note  that  once  an 
inclination  change  is  initiated,  EVA  maintains  a  constant  a  magnitude  throughout  the  transfer.  The 
modelling  equations  assume  that  the  direction  of  a  is  reversed  at  the  maximum  and  minimum  orbital 
latitude  positions  (approximately  90°  away  from  the  nodes)  to  maintain  inclination  change  in  the  same 
direction.  If  the  desired  inclination  change  is  not  achieved  at  the  end  of  the  transfer,  the  entire  transfer  is 
iterated  using  a  secant  method  to  converge  on  the  appropriate  value  of  a  in  order  to  meet  the  desired  final 
inclination. 
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The  code  permits  delaying  out-of-plane  thrusting  until  the  vehicle  reaches  a  specified  altitude.  This 
provides  the  option  to  “quickly”  raise  the  orbital  altitude  beyond  the  majority  of  atmospheric  interference 
and  lengthy  eclipses  before  diverting  available  thrust  for  an  inclination  change. 

Atmospheric  Drag.  The  drag  force  acting  on  the  vehicle  is 
computed  using  the  familiar  equation 

D  =  Cdq  A  (24) 

where  Cd  is  the  drag  coefficient,  A  is  the  reference  cross  sectional 
area  of  the  vehicle,  and  where  q  is  the  dynamic  pressure  given  by 

q  =  ^P  V2  (25) 

During  an  actual  mission,  the  vehicle’s  solar  arrays  will  be 
continuously  tracking  the  sun  and  thus  not  be  pointing  at  a  constant 
angle  with  respect  to  the  vehicle’s  velocity  vector.  This  means  that  the 
vehicle’s  reference  cross  sectional  area  is  constantly  changing.  As  the  arrays  move  from  90°  to  0°  with 
respect  to  the  velocity  vector,  the  reference  area  changes  from  100%  of  the  total  cross  sectional  area  to  0%. 
EVA  calculates  the  worst  case  scenario  where  the  reference  area  equals  100%  of  the  total  cross  sectional 
area  at  all  times.  Circular  orbital  velocity  at  the  time  of  the  drag  calculation  is  used  to  compute  dynamic 
pressure.  A  self-contained  subroutine,  incorporating  the  1976  U.S.  Standard  Atmosphere,  is  used  to 
compute  atmospheric  density  as  a  function  of  altitude.8  This  is  a  static  atmosphere  model.  The  drag 
coefficient  is  available  as  an  input  to  the  program.  The  drag  force  is  then  factored  into  Equation  6, 
effectively  diminishing  the  available  in-plane  thrust 

The  effects  of  atmospheric  drag  on  the  trajectory  is  assumed  negligible  above  1,000  kilometers 
altitude,  and  is  not  considered  once  the  trajectory  reaches  that  altitude.  At  the  lower  altitudes,  the 
assumption  is  made  that  the  vehicle’s  available  in-plane  thrust  to  drag  ratio  must  exceed  1.1  in  order  to 
continue  on  the  trajectory.  That  constraint  is  imposed  by  computing  an  effective  thrust,  Teff,  which  is 
constrained  to  be  greater  than  zero. 

Teff=  0.9  cosa  F  -  D/fs  (26) 

where  fs  represents  the  fraction  of  the  orbit  the  vehicle  is  in  sunlight  and  capable  of  providing  thrust. 

Power  System  Interaction 

Occultation.  To  determine  the  amount  of  time  during  which  the  thrusters  will  not  receive  power 
from  the  power  system,  the  fraction  of  the  transfer  vehicle’s  orbit  spent  in  the  earth’s  occultation  zone 
must  be  calculated.  The  eclipse  duration  is  a  function  of  the  date,  which  determines  the  earth-sun 
geometry,  as  well  as  the  orbital  parameters  of  the  spacecraft.  From  spherical  trigonometry  and  Figure  8, 

sin  5  _  sin[(os  (t  -  Q] 

sin  e  sin  90°  nn\ 
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Autumn 


Summer 


and 


8  =  sin  sin  e  sin[(os(t  - 1^]  j 

cosjw  ,|t  -  t0jj  =  cosa  cos8  +  sina  sin5  cos90° 

a  =  cos'1!— —^t— -I 
(  cos  5  | 


(28) 


(29) 


Knowing  a  and  5  and  the  spacecraft’s  orbital  parameters,  the  angle  between  the  sun  and  the  orbital 
plane  can  be  calculated.9 


P  =  sin'^cos  8  sin  i  sin(f2  -  a)  +  sin  8  cos  i} 


(30) 


The  spacecraft’s  right  ascension  of  ascending  node  will  slowly  rotate  during  the  transfer  due  to  the 
earth’s  oblateness.  This  rotation  will  affect  the  beta  angle;  therefore,  the  eclipse  duration.  The  regression 
of  the  node  can  be  calculated  using  Equation  31. 

Q  =  J2(y^cosi 

2\rl  (31) 
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Using  Figure  3  and  defining  the  angle  (5*,  the  fraction  of  the  orbit  spent  in  eclipse  can  be  defined  as  in 
Equation  29.  Note  that  the  fraction  of  the  orbit  spent  in  sunlight  (fs)  is  simply  1-  fe.  This  fraction  can  then 
be  incorporated  into  Equation  8  to  calculate  the  orbit  transfer  possible  with  the  given  geometry. 


Figure  9 

Sun-Earth  Occultation  Geometry 


Solar  Array  Power  Degradation.  The  slow  spiral  transfer  through  the  Van  Allen  radiation  belts 
can  result  in  significantly  degraded  solar  array  power  output  and  reduced  vehicle  performance.  In 
conjunction  with  the  orbit  transfer  calculations,  the  program  automatically  computes  the  cumulative  fluence 
experienced  by  the  vehicle  and  determines  the  resulting  power  degradation  for  the  specified  solar  cell  type. 
Logic  within  the  program  permits  the  degradation  calculations  to  be  bypassed  as  well  as  automatically 
terminating  the  run  should  power  degrade  to  zero. 

For  each  transfer  orbit  revolution,  the  annual  equivalent  1  MeV  electron  fluences  from  trapped 
electrons  and  protons  are  determined  from  tabular  fluence  data  stored  in  the  program.10  The  data  is 
tabulated  for  varying  orbital  altitudes,  inclinations  (up  to  30°),  and  shield  thicknesses  (up  to  60  mils). 
Linear  interpolation  is  used  to  compute  the  fluence  levels  for  specified  shield  thicknesses  at  the  current 
transfer  orbit  altitude  and  inclination.  The  combined  front  and  back  side  fluence  experienced  by  the  array 
cells  is  computed  using  Equation  33.  Different  shield  thicknesses  can  be  specified  for  front  and  back 
sides. 

^combined  =  [^electron  +  0proton]  "*■  [^electron  +  9proton)]  (33) 

Front  Side  of  Array  Back  Side  of  Array 

The  combined  fluence  (for  the  front  and  back  sides  of  the  solar  array)  is  numerically  integrated  to 
yield  the  cumulative  fluence  which  the  vehicle  arrays  experience  during  the  transfer.  This  cumulative 
fluence  value  is  then  translated  into  percent  power  degradation  using  tabulated  data  and  logarithmic 
interpolation.  The  degraded  power  value  is  used  in  equations  17  and  18  to  lower  thrust  and  delivered  Isp 
as  appropriate  for  the  specific  electric  propulsion  system  design. 
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System  Sizing.  In  addition  to  the  orbital  mechanics  calculations,  a  mission  planning  tool  needs 
calculations  to  determine  weights  of  payload  and  necessary  subsystems.  The  system  sizing  logic  in  EVA 
takes  the  actual  EOTV  dry  mass  and  divides  it  into  five  major  categories:  arrays,  structure,  thermal 
management,  propulsion  system,  and  payload.  The  user  is  allowed  to  input  specific  mass  (mass  per  unit 
power)  relationships  for  the  five  categories  in  order  to  determine  system  weight.  Also  input  are  specific 
power  (power  per  unit  mass)  and  total  system  power  level,  which  are  used  to  calculate  the  necessary  solar 
array  surface  area  mandatory  for  drag  calculations.  Many  different  array  types  can  be  modelled  by 
changing  the  specific  power  parameter. 


Figure  10 
EOTV  Subsystems 

The  propulsion  system  and  power  system  masses  are  determined  simply  by  multiplying  their 
respective  specific  masses  and  the  system  power. 


The  thermal  management  system  mass  is  dependent  on  the  power  processing  unit  (PPU)  efficiency. 
The  thermal  management  system  must  reject  the  power  which  the  PPU  cannot  convert  into  usable  thrust. 
Heat  generated  by  thruster  inefficiency  is  considered  to  be  handled  by  direct  radiation  from  the  thruster  to 
space  and,  hence,  does  not  affect  the  mass  of  the  thermal  rejection  hardware.  The  thermal  management 
system  mass  (Mt)  is  determined  by  the  following  equation 

M.=(1-v)*p>*w.  (34) 

where  Pi  is  the  initial  system  power  and  Wt  is  the  thermal  management  system  specific  mass. 
Structure  mass  is  calculated  using  a  fixed  percent  of  the  dry  mass  of  the  EOTV  and  its  payload.  This 
percentage  is  a  user-defined  entry  in  the  input  file — currently  the  default  value  is  12%  which,  after  talking 
with  several  people  in  our  in-house  structure  group,  was  decided  upon  as  a  good  starting  estimate. 

The  payload  is  everything  not  already  accounted  for  in  the  above  categories.  It  is  determined  by 
subtracting  the  sum  of  the  above  categories  from  the  total  EOTV  dry  mass. 
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Typical  Problems 

In  its  current  form  EVA  can  provide  mission  planners  with  useful,  reasonably  accurate  answers  to  a 
wide  variety  of  questions.  Even  without  internal  logic  for  paramaterization  it  is  simple  to  generate  a  quick 
look  of  parametric  effects  on  any  given  system.  One  important  figure  of  merit,  for  example,  is  trip  time 
for  delivery  of  a  satellite.  EVA’s  flexibility  allows  the  user  to  determine  effects  on  trip  time  due  to  power 
levels,  array  protection,  propulsion  system  efficiency,  and  drag,  among  others.  It  will  also  allow  for  the 
determination  of  optimal  deployment  altitudes-trading  the  increased  throw  weight  by  the  launch  vehicle  to 
lower  altitudes  with  the  lower  accelerations  possible  at  these  lower  altitudes  due  to  the  more  severe  drag 
environment.  The  cumulative  radiation  fluence  calculation  also  provides  an  opportunity  to  study  the 
effects  of  radiation  on  the  payload.  The  sizing  logic  allows  the  user  to  vary  subsystem  specific  masses  for 
subsystems  such  as  arrays,  power  processors,  thermal  management,  and  structure  and  to  understand  the 
payload  capability  for  any  number  of  EOTV  designs.  This  flexibility,  coupled  with  the  program’s 
computational  efficiency,  make  it  an  extremely  useful  preliminary  mission  analysis  tool,  capable  of  quickly 
studying  a  variety  of  missions.  One  such  mission  is  shown  in  the  following  example  problem. 

A  current,  chemical  OTV  that  is  used  on  the  Titan  IV  is  the  Inertial  Upper  Stage  (IUS).  It  can  deliver 
2,268  kg  to  geosynchronous  orbit  using  12,429  kg  of  propellant  and  with  a  dry  weight  of  2,279  kg.  How 
much  of  the  propellant  weight  could  be  saved  by  using  an  EOTV  to  deliver  the  same  2,268  kg  to  the  same 
geosynchronous  orbit? 

To  use  EVA  to  answer  this,  first  go  to  the  EVA  input  file.  Check  the  initial  and  final  orbit  altitudes  to 
make  sure  they  reflect  a  LEO  to  GEO  transfer.  Also  check  to  see  that  the  inclination  change  and  the 
altitude  to  start  changing  inclination  are  correct.  The  drag  coefficient  is  defaulted  to  2.0;  if  more  accurate 
data  is  available,  then  this  variable  may  be  input  directly.  The  day  of  the  year  desired  to  launch  is  also  an 
input  since  it  affects  the  amount  of  time  spent  in  the  earth’s  occultation  zone.  The  default  value  is  80, 
corresponding  to  a  vernal  equinox  launch,  and  need  not  be  input  unless  a  specific  launch  day  is  desired. 
Check  the  flags  for  power  degradation  and  atmospheric  drag.  Enter  the  data  for  the  specific  type  of  EOTV 
you  wish  to  run.  For  purposes  of  this  example  the  values  used  are  typical  of  an  arcjet  type  electric  thruster 
and  are  shown  in  Table  1.  Appendix  A  contains  the  complete  input  variable  list  with  corresponding 
definitions. 


Table  1.  Example  Problem  Inputs 


Isp . 1000  sec 

system  power . 30  kW 

vehicle  cross  sectional  area 

(excluding  solar  arrays) . ...5  sq  m 

PPU  efficiency . 0.9 

thruster  efficiency . 0.38 

initial  fueled  vehicle  mass . 5,541  kg 

specific  mass  of  propulsion  system . 3.2  kg/kw 

specific  mass  of  thermal  management  system. 18  kg/kw 
structure  fraction . 0.14 _ 
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Now  close  the  input  file  and  run  the  program.  After  several  seconds,  the  output  file  will  be  completed 
(see  Appendix  A  for  the  full  output  listing). 

Table  2.  Partial  Output  Listing 


ELECTRIC  PROPULSION  MISSION  SUMMARY 

Structural  mass:  431.20 

Thermal  management  mass:  54.00 
Propulsion  system  mass:  96.00 

Power  mass:  230.77 

811.97  *Total  EOTV  dry  mass 

Payload:  2268.03 

3080.00  **Total  Ma33  0  GEO 

Propellant:  2460.77 

5540.77  ***Total  Mass  0  LEO 

MISSION: 

Initial  Altitude  (km):  300. 

Final  altitude  (km)  :  35775. 

Required  Delta-V  (m/s):  6104.76 

Time  in  sun  (days):  183.75 

Time  in  shadow  (days):  36.94 

Total  trip  time(days):  220.69 

Total  number  of  orbits:  1522.27 


In  answer  to  the  problem,  from  the  output  we  see  that  this  EOTV  would  weigh  approximately  812  kg, 
use  2,460  kg  propellant,  and  carry  2,268  kg  payload  to  geosynchronous  orbit.  The  trip  would  take  221 
days  of  which  37  days  are  spent  in  the  earth’s  shadow  (non-thrusting).  This  means  that  9,970  kg  of 
propellant  could  be  saved  by  using  an  EOTV.  This  is  important  because  the  weight  savings  is  enough  to 
be  able  to  launch  the  vehicle  from  an  Atlas  II  rather  than  a  Titan  IV  resulting  in  a  potential  savings  of  many 
millions  of  dollars. 

A  useful  feature  of  the  code  is  its  ability  to  recognize  the  fact  that  atmospheric  drag  may  be  greater 
than  the  thrust  of  the  vehicle.  In  such  cases,  the  code  will  automatically  raise  the  starting  altitude  up  to 
such  a  height  that  thrust  exceeds  drag  by  at  least  10%.  Similarly,  if  during  the  vehicle’s  trek  through  the 
Van  Allen  radiation  belts,  the  solar  arrays  become  so  damaged  by  radiation  that  the  power  produced  falls  to 
zero,  the  code  will  alert  the  user  to  that  fact  and  terminate  the  run  (Appendix  A). 

In  The  Future... 

As  previously  mentioned,  one  of  the  objectives  in  creating  EVA  was  to  make  it  easy  to  update  to 
include  better  models,  and  more  user-friendly  features.  There  are  several  upgrades  on  the  horizon  for 
future  versions  of  the  code.  Two  of  these  proposed  features  are  automatic  parametric  runs  and 
optimization. 
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Parametric  studies  are  often  essential  in  preliminary  mission  planning  since  many  of  the  vehicle  and 
mission  variables  are  not  defined.  A  useful  feature  of  the  EVA  would  be  the  capability  to  run  a  range  of 
values  of  specific  variables  in  a  single  run,  and  obtain  a  graph  of  the  results  (e.g.  trip  time  vs  system 
power).  Of  course,  parametric  runs  are  currently  possible  with  EVA  by  simply  editing  the  input  file  and 
changing  the  appropriate  variable(s)  and  rerunning  the  code.  This,  however,  is  not  as  convenient  as 
having  the  code  do  that  automatically. 

Optimization  is  another  userful  feature.  Many  of  the  variables  directly  impact  other  variables.  For 
instance,  more  system  power  is  desirable  because  it  increases  thrust,  thus  reducing  trip  time,  but  it  also 
increases  weight  and  drag  which  increases  trip  time.  Within  a  certain  range  there  is  an  optimum  power 
level  which  minimizes  trip  time.  An  optimization  routine  would  be  useful  in  finding  such  optimal  values. 

SUMMARY/CONCLUSION 

Recent  activity  in  mission  planning  for  EOTV  systems  and  flight  experiments  such  as  ELITE,  led  to 
the  need  for  a  mission  planning  tool  which  incorporated  several  performance  and  system  level  parameters. 
Existing  tools  within  government  agencies  did  not  adequately  fulfill  that  need  and,  hence,  the  Electric 
Vehicle  Analyzer  was  created.  Although  EVA  does  not  model  detailed  attitude  and  flight  control 
algorithms  it  will  perform  EOTV  mission  studies  with  enough  accuracy  to  impact  vehicle  and  mission 
design  while  running  efficiently  and  quickly  in  a  VAX  or  PC  environment.  The  code  allows  the  user  a 
significant  degree  of  flexibility  in  input  parameters,  allowing  the  capability  to  analyze  many  different 
missions  performed  by  many  different  vehicles. 

EVA  has  become  a  valuable  part  of  the  mission  planning  toolbox  at  the  Astronautics  Laboratory.  It 
will  continue  to  be  updated  and  refined  while  always  keeping  an  eye  toward  efficiency  and  user  flexibility. 
The  growing  national  interest  in  fielding  EOTV  systems  ensure  that  the  program’s  capabilities  will  be  well 
exercised  in  the  years  to  come.  The  development  of  the  code  focused  on  maintaining  the  flexibility  to 
consider  a  number  of  types  of  missions  and  vehicles.  The  authors  welcome  any  comments  from  the 
electric  propulsion  community  as  to  things  we  may  have  missed  or  suggestions  as  to  how  the  program  can 
be  made  more  responsive  for  specific  applications. 
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APPENDIX  A 


INPUT  FORMAT 

The  input  variables  are  entered  into  the  program  in  the  form  of  a  namelist  input  file.  This  method  allows 
the  user  to  change  any  or  all  of  these  variables  between  runs  in  a  timely  fashion  since  changing  the  input 
file  does  not  require  any  recompilation.  Table  3  provides  a  definition  of  each  input  variable  and  its  default 
value  if  any. 


Table  3.  Definition  of  Program  Input  Variables 


ALPHA 

Guess  for  out-of-plane  bum  angle  used  to  change  inclination  (will  be  adjusted  for 
INCLF),  degrees 

ALTI 

Initial  orbit  altitude,  km 

ALTF 

Final  orbit  altitude,  km 

ALTINCL 

Altitude  above  which  orbit  inclination  change  is  performed,  km  (ALTINCL  must 
be  set  below  ALTF  if  there  is  an  inclination  change.) 

CD 

Coefficient  of  Drag  (Default  value  =  2.0) 

DATE 

Date  during  the  year  when  the  transfer  is  initiated  in  days  from  year  start.  (Default 
value  =  80.0) 

ENGEFF 

Thruster  efficiency 

FAREA 

Factor  for  adjusting  array  cross-sectional  area  in  drag  calculations  to  account  for 
varying  orientation  with  respect  to  the  velocity  vector.  (Default  =  1.00) 

FDENS 

Multiplication  factor  for  atmospheric  density  term  in  drag  computations.  Used  to 
approximate  worst  case  drag  conditions.  (Default  value  =  1.00) 

ICELL 

Index  specifying  solar  cell  type  for  determining  array  power  degradation  (Input 
for  EDEGRD  =  1) 

1  -  Si  BBSFR  Thin  Cell,  Gridded  Back,  3.4  mil 

2  -  ASEC  OMCVD  GaAs/Ge 

IDRAG 

Flag  for  including  effect  of  atmospheric  drag 

0  -  No  drag  effects  included 

1  -  Drag  effects  included  using  1976  Std.  Atmosphere 

EDEGRD 

Flag  ior  including  solar  array  power  degradation  effects 

0  -  Neither  cum  fluence  or  array  power  degradation  is  computed  during 
transfer 

1  -  Cum  fluence  computed  but  array  power  not  degraded;  power  output  and 
thrust  remain  constant 

2  -  Cum  fluence  computed  and  array  power  degraded 

IPLOT 

Rag  for  creating  plot  output  File  during  run  (Default  =  0) 

0  -  No  plot  output  file  created 

1  -  Create  plot  file  (Unit  =  8)  during  run 
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IPRINT 


INCLI 


INCLF 


ISP 


NPRINT 


NPLOT 


OMEGAI 


PROPSYS 


PPUEFF 


PWRDEN 


PWRWGT 


SCAREA 


SCMASS 


STRFRACT 


SYSPWR 


THMGMT 


THSTPWR 


TSHIELD(l) 


TSHIELD(2) 


of  output  to  be  generated  during  run 


r  data  only  (Default) 


0-S 


1  -  Sum 


Initial  orbit  inclination,  degrees  (Default  =  28.5  deg) 


Final  desired  orbit  inclination,  degrees 


Thruster  specific  impulse,  sec 


Orbit  interval  at  which  trajectory  data  will  be  output  if  IPRINT  >  0. 
(Default  value  =  100) 


Orbit  interval  at  which  trajectory  data  will  be  written  to  plot  tape  if  IPLOT  >  0. 
(Default  value  =  100) 


Right  Ascension  of  Ascending  Node  at  start  of  transfer,  degrees.  (Value  will  be 
adjusted  during  transfer  to  account  for  perturbation  effects  of  equatorial  bulge.) 


Specific  mass  of  propulsion  system,  k 


Efficiency  of  the  power 


Power  density  of  solar  artays,  W/m2 


Specific  power  of  solar  arrays,  W/k 


Spacecraft  body  cross-section  excluding  solar  array  area,  m2. 


Initial  fueled  mass  of  spacecraft  excluding  solar  arrays,  k 


Structural  fraction  of  vehicle 


Initial  system  power  provided  by  solar  arrays,  kW.(Remains  constant  if  IDEGRD 
=  0)  provided  by  solar  arrays,  kW 


Specific  mass  of  thermal  management  unit,  kg/kW 


Power  required  per  electric  thruster,  kW 


Array  front-side  effective  shield  thickness  for  use  in  determining  fluences,  0-60 
mils  (input  for  IDEGRD  =  1) 


Array  back-side  effective  shield  thickness  for  use  in  determining  fluences,  0-60 
mils  (input  for  IDEGRD  =  1) 


OUTPUT  FORMAT 

The  example  problem  output  is  shown  below.  The  names  of  the  variables  output  are  printed  at  the  top  of 
the  page  as  a  header,  and  then  the  values  of  these  variables  are  printed  as  a  block  at  successive  time 
intervals.  At  the  end  of  this  trajectory  output,  a  summary  of  weights  and  propulsion  requirements  is 
printed. 


Table  4.  Example  Problem  Output 

ELECTRIC  VEHICLE  ANALYZER  (EVA)  PROGRAM 


TIKE  (DAYS) 

ALT  (ICM) 

INCL  (DEC) 
KUM  ORBITS 

VEL  (M/JEC) 
ROOT  (M/SEC) 
GAtHA  (DEG) 

TOT  DV  (M/S) 
PERIOD  (HRS) 
BETA  (DEG) 

CUM  FLUES  IMEV) 
FLUES  (MEV/TR) 
ECLIPSE  FRACT 

POWER  DEGRAD 
POWER  (KM) 
ATM  DRAG  (N) 

TCTT  THRUST  (N) 
ALPHA  (DEG) 

P  LA  THRUST  (M) 

S/C  MASS  (KCI 
PROP  MASS  (KG) 
SPEC  -MP  (SEC) 

6.395 

463.603 

28.0919 

100.00 

7632.83 

0.3379 

0.00254 

133.53 

1.5643 

-19.1834 

0.1031D+11 

0.1043D+13 

0.37509 

1.00000 

30.0000 

0 . 1339D-01 

2.092555 

40.0362 

1.602139 

5238.1854 

71.8146 

1000.00 

13.060 

671.077 

27.6319 

200.00 

7519.67 

0.3790 

0.00289 

281.94 

1.6359 

-22.0584 

0.8070D+11 

0 .84930+13 
0.34828 

1.00000 

30.0000 

0.5978D-03 

2.092555 

40.0362 

1.602139 

5159.5026 

150.4974 

1000.00 

20.048 

909.693 

27.1210 

300.00 

7395.54 

0.4114 

0.00319 

444.16 

1.7196 

-6.6983 

0.3824D+1 2 
0.2868D+14 
0.33825 

1.00000 

30.0000 

0. 7143D-04 

2.092555 

40.0362 

1.602139 

5074.8578 

235.1422 

1000.00 

27.414 

1185.871 

26.5505 

400.00 

7259.27 

0.4594 

0.00363 

622.21 

1.818'’ 

15.4963 

C.  ..4.30+13 

0. 1002D+15 
0.31184 

0.99718 

29.9155 

0.00000 

2.086658 

40.0362 

1.597791 

4983.5456 

326.4544 

1000.00 

35.237 

1521.700 

25.8842 

500.00 

7103.30 

0.5385 

0.00435 

826.01 

1.9405 

34.1077 

O.SS'.iO*  *> 

0.3166D+15 

0.24780 

0.98119 

29.4357 

0.00000 

2.053197 

40.0362 

1.572359 

4881.0393 

428.9607 

1000.00 

43.643 

1948.813 

25.0763 

600.00 

6918.73 

0.6305 

0.00523 

1067.21 

2.0998 

40.5074 

0.2024D+14 

0.1104D+16 

0.17969 

0.94982 

28.4946 

0.00000 

1.987552 

40.0362 

1.522378 

4762.4487 

547.5513 

1000.00 

52.797 

2465.812 

24.1520 

700.00 

6713.45 

0.6675 

0.00570 

1335.47 

2.2983 

33.5985 

0. 7136D+14 
0.3344D+16 
0.18770 

0.90273 

27.0820 

0.00000 

1.889020 

40.0362 

1.447111 

4633.9303 

676.0697 

1000.00 

62.851 

3053.196 

23.1652 

800.00 

6501.04 

0.6846 

0.00604 

1613.08 

2.5309 

21.5821 

0.2255D+15 

0 . 8290D+1 6 
0.20917 

0.83926 

25.1779 

0.00000 

1.756207 

40.0362 

1.345664 

4504.5838 

805.4162 

1000.00 

73.957 

3724.050 

22.1106 

900.00 

6281.47 

0.7153 

0.00653 

1900.05 

2.8056 

10.9755 

0. 5810D+15 

0 . 1503D+1 7 
0.21031 

0.76932 

23.0796 

0.00000 

1.609848 

40.0362 

1.233623 

4374.6742 

935.3258 

1000.00 

86.329 

4518.829 

20.9481 

1000.00 

6048.06 

0.7754 

0.00735 

2205.12 

3.1428 

4.7472 

0.1199D+16 

0.2045D+17 

0.19769 

0.71050 

21.3151 

0.00000 

1.-186771 

40.0362 

1.139169 

4240.6710 

1069.3290 

1000.00 

100.304 

5513.106 

19.6075 

5789.69 

0.8772 

2542.86 

3.5822 

0.2008D+16 

0 . 2067D+1 7 

0.66634 

19.9903 

1.394358 

40.0362 

4097.1034 

1212.8966 
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1100.00 

0.00869 

3.5645 

0.17943 

0.00000 

1.068151 

1000.00 

116.452 

6852.723 

17.9678 

1200.00 

5488.77 

1.0536 

0.01102 

2936.29 

4.2035 

6.1827 

0.2796D+16 

0 . 1447D+17 
0.15691 

0.63803 

19.1408 

0.00000 

1.335106 

40.0362 

1.022549 

3935.9795 

1374.0205 

1000.00 

135.878 

8881.845 

15.7744 

1300.00 

5110.84 

1.3876 

0.01560 

3430.57 

5.2047 

9.8067 

0.3311D+16 

C.5507D+16 

0.12688 

0.62356 

18.7067 

0.00000 

1.304829 

40.0362 

0.999162 

3742.5053 

1567.4947 

1000.00 

161.259 

12675.828 

12.3550 

1400.00 

4573.79 

2.1385 

0.02691 

4133.45 

7.2550 

9.8986 

0.3466D+16 

0.4148D+15 

0.09456 

0.61966 

18.5897 

0.00000 

1.296665 

40.0362 

0.992790 

3483.6400 

1826.3600 

1000.00 

203.230 

25005.406 

4.6179 

1500.00 

3563.84 

5.2984 

0.08638 

5458.91 

15.2626 

-4.0065 

0.3488D+16 

0 . 1629D+15 
0.06193 

0.61911 

18.5734 

0.00000 

1.295525 

40.0362 

0.991914 

3043.2086 

2266.7914 

1000.00 

220.694 

35775.000 

-0.0304 

1522.27 

3075.06 

9.4586 

0.17761 

6104.76 

6.4168 

-13.8908 

0.3492D+16 

0.2532D+14 

0.00000 

0.61901 

18.5703 

0.00000 

1.295311 

40.0362 

0.991740 

2849.2350 

2460.7650 

1000.00 

ELECTRIC  PROPULSION  MISSION  SUMMARY 


THRUSTERS 


Thrust/Thruster  (N) :  0.43 
Quantity:  3 
Total  Thrust  (N) :  1.30 
Specific  Impulse  (sec):  1000. 


POWER: 

Power  Density  (W/M**2) :  130.00 
Specific  Power  (W/kg) :  130.00 

Total  Power  (kW) :  30.00 

Excess  Power  (kW) :  0.00 

Array  Area  (M**2) :  230.77 

Array  Mass  (kg):  230.77 

System  Efficiency:  0.342 

SPACECRAFT  WEIGHT  BREAKDOWN  (kg) : 

Structural  mass:  431.20 

Thermal  management  mass:  54.00 


Propulsion  system  mass:  96.00 

Power  mass:  230.77 

811.97 

♦Total  EOTV  dry  mass 

Payload:  2268.03 

3080.00 

♦♦Total  Ma3s  @  GEO 

Propellant:  2460.77 

5540.77 

♦♦♦Total  Mass  @  LEO 
MISSION: 


Initial  Altitude  (km) :  300 
Final  altitude  (km) :  35775 
Required  Delta-V  (m/s):  6104.76 

Time  in  sun  (days):  183.75 
Time  in  shadow  (days):  36.94 
Total  trip  time(days):  220.69 
Total  number  of  orbits:  1522.27 


A  useful  feature  of  the  code  is  its  ability  to  recognize  the  fact  that  atmospheric  drag  may  be  greater 
than  the  thrust  of  the  vehicle.  In  such  cases,  the  code  will  automatically  raise  the  starting  altitude  To  such  a 
height  that  thrust  exceeds  drag  by  at  least  10%. 


Similarly,  if  during  the  vehicle’s  trek  through  the  Van  Allen  radiation  belts,  the  solar  arrays  become 
so  damaged  by  radiation  that  the  power  produced  falls  to  zero,  the  code  will  alert  the  user  to  that  fact  and 
terminate  the  run.  Table  4  demonstrates  this  type  of  run. 
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Table  5.  Low  Altitude/High  Power  Degradation  Case 


*  DRAG  EXCEEDED  90*  OF  IN-PLANE  THRUST  FOR  INPUT  PARAMETERS;  STARTING  ALTITUDE  RAISED  FROM  160.900  KM  TO  260.900  KM 


ELECTRIC  VEHICLE  ANALYSIS  (EVA)  PROGRAM 


TIME  (DAYS)  ALT 

(KM) 

VEL 

(M/SEC) 

TOT  DV 

(M/S) 

CUM  FLUEN  (MEV) 

POWER  DEGRAD 

TOT 

THRUST 

(N) 

S/C  MASS  (KG) 

INCL 

(DEG) 

RDOT 

(M/SEC) 

PERIOD 

(HRS) 

FLUEN 

(MEV/YR) 

POWER  (KW) 

ALPHA  (DEG) 

PROP 

MASS (KG) 

NUM 

ORBITS 

GAMMA 

(DEG) 

BETA 

(DEG) 

ECLIPSE 

FRACT 

ATM  DRAG  (N) 

PLN 

THRUST 

(N) 

SPEC 

IMP (SEC) 

3.123 

284.011 

7735.03 

24.44 

-0.6639D+04 

1.00000 

2.007303 

13084.4235 

0.0000 

0.1042 

1.5038 

0.9641D+06 

28.0000 

45.0000 

32.6565 

50.00 

0.00077 

1.1996 

0.40679 

0.3988D+00 

2.007303 

1000.00 

6.267 

316.646 

7716.15 

49.30 

0.2133D+05 

1.00000 

2.007303 

13051.2999 

0.0000 

0.1333 

1.5147 

0.5628D+07 

28.0000 

45.0000 

65.7801 

• 

100.00 

0.00099 

2.4288 

0.40173 

0.1984D+00 

2.007303 

1000.00 

888.759 

11805.240 

4682.00 

3089.34 

0.3353D+19 

0.03145 

0.063123 

9572.3149 

0.0000 

0.0460 

6.7777 

0.1758D+19 

0.8805 

45.0000 

3544.7651 

5600.00 

0.00056 

9.4204 

0.10184 

0 .0000D+00 

0.063193 

1000.00 

902.912 

11859.786 

4675.00 

3096.35 

0.3420D+19 

0.02974 

0.059691 

9565.4774 

0.0000 

0.0433 

6.8082 

0.1735D+19 

0.8326 

45.0000 

3551.6026 

5650.00 

0.00053 

4.0971 

0.11152 

O.OOOOD+OO 

0.059759 

1000.00 

914.013 

11900.400 

4669.80 

3101.55 

0.3473D+19 

0.02843  0.057077 

9560.4094 

0.0000 

0.0414 

0.4442 

0. 171 7D+1 9 

0.7962  45.0000 

3556.6706 

5689.07 

0.00051 

-0.3410 

0.11344 

0 .0000D+00  0.057081 

1000.00 

•  ZERO 

POWER  AVAILABLE 

FOR  THRUSTERS  AT  11900.446  KM  DUE 

TO  ARRAY  DEGRADATION 

;  RUN  TERMINATED 
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APPENDIX  B 


PROGRAM  EVA 

IMPLICIT  DOUBLE  PRECISION  (A-M.O-Z) 

CHARACTER  WEPTYP*(18),WPWRTYP*(50) 

INTEGER  US,roRAG,IDEGRD,ICELLJST(2),IPL0TJPRINT 
INTEGER  PWRTYP.EPTYP 

DIMENSION  FLUEN(5),TSHIELD(2),TSHLDR(8),FST(2) 

COMMON/CONST/MU, SQMUJRE,PI,PI2,DTS,DTR,KTM,KTM3,G0,r YEAR.SEPSLN 
COMMON/FLAGS/IDRAG,IDEGRD4CELL,IST  JFST  JPRENT  JPLOT  .NPRINT 
COMMON/TRAJ/T,DAYSJDATEALTtALTI,ALTF,VELtRDOT,INCL,INCLUNCLF, 

*  ALTINCL,CALPH,OMEGAI,BETA,FE,CD,AREA,FDENS,GAMMA,ORBITS, PERIOD, 

*  FPWRJFLUEN.CFLUEN, THRUST  ,THRUSTI,THRUSTP,DRAG,ISP,ETA,DVTOT, 

*  DRGLOSS,GRVLOSS,MASSI,MASS,MPROP,SYSPWR£PPWRJPOWER,TECLPS( 

*  NUMTHST 

COMMON/MASS/STR,THMGMT,INEFF,PROPSYS,EPPWERlSASM,MEOTV,MPAY,AAREA, 

*  AMASS,XSPWRJPWRWGT,PWRDEN,SCMASS,PPUEFF,ENGEFF,STRFRACT 
COMMON/SUNyTSUN 

Initialize  Program  Constants  and  Variables 

DATA  MU.SQMU.RE/3.98600800D14,  1.9964989D7, 6.37813500D6/ 

DATA  PI J*I2 JDTS/3 .141 5926535898 DO,  6.2831 85307 1796D0,  8.64D4/ 

DATA  DTR.KTM.KTM3/0. 174532925 199433D-1,  1000.D0,  1.D9/ 

DATA  G0JDYEAR£PSILON/9.806194D0,  365.25D0,  23.4432D0/ 

DATA  IDRAG JDEGRD.ICELL JST.IPLOT JPRINT/1 ,1 , 1 ,2*  1 ,0,0/ 

DATA  DATE ,FDENS,CD,INCLI,INCLF/80.D0,1  .DO,2.0DO,2*28.5DO/ 

DATA  TSHLDR/0.DO,l.DO,3.DO,6.DO,12.DO,2O.DO,3O.DO,6O.DO/ 

DATA  TSHIELD  JLUEN  JST^ARE A/2*O.DO^*O.DO,2*O.DO,1  .ODO/ 

namelist/evain/altialtfjnclunclf.altinclalpha.omegai.date. 

*  IDRAG,SCAREA,CD,FDENSJrAREA,IDEGRD,TSHIELD,ICELL.SCMASS.SYSPWR, 

*  ISP,ENGEFF,PPUEFF,NTHSTR,THSTPWR,PWRDEN,PWRWGT,IPRINT,NPRINT, 

*  IPLOT ,NPLOT,THMGMT,PROPS  Y S  .STRFR ACT 
OPEN  (UNIT=5,  STATUS='OLD',  FELE=  INPUT) 

OPEN  (UNIT=6,  STATUS^’NEW’,  FELE=E  V  AOUT) 

IF(IPLOT.EQ.l)  OPEN  (UNIT=8,  STATUS=’NEW,  FTLE='EV  APLOT) 

Read  Inputs  and  Perform  Initial  Computations 

READ(5,EVAIN) 

WRITE(6  JEV  AIN) 

INCLI  =  DTR*INCLI 
INCLF  =  DTR*INCLF 
CALPH  =  DCOS(DTR*  ALPHA) 

SEPSLN  =  DSIN(DTR*EPSILON) 

ETA  =  ENGEFPPPUEFF 


Determine  No.  of  Thrusters  and  Total  System  Thrust  at  Mission  Start 

NUMTHST  =  DINT(SYSPWR/THSTPWR) 

EPPWR  =  FLOAT(NUMTHST)*THSTPWR 
XSPWR  =  SYSPWR  -  EPPWR 
MDOT  =  (2.0D3*ETA*EPPWR)/(GO*ISP)**2.0D0 
THRUSH  a  GO*ISP*MDOT 


Determine  S/C  area  and  mass  including  arrays... 
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ARRAY=1 
AMASS  =  0.D0 

S  ASM=(1/PWRWGT)*  1000.0D0 
MASSI  =  AMASS  +  SCMASS 
AAREA  =  S  YSPWR*  lOOO.DO/PWRDEN 
AREA  =  AAREA  +  SCAREA 

Determine  Fluence  Interpolation  Factors  for  Specified  Front  and 
Back-Side  Array  Shield  Thicknesses  (60  mils  Maximum  Thickness) 

DO  30  IS  =  U 
DO  10 1  =  2,8 

IF  (TSHIELD(IS)l-T.TSHLDR(I))  GOTO  20 
10  CONTINUE 
IST(IS)  =  8 
FST(IS)  »  0.D0 
GOTO  30 

20  ISTaS)  =  I-l 

FST(IS)  =  (TSHIELD(IS)  -  TSHLDR(I- 1))/(TSHLDR(I)  -  TSHLDR(1-1» 

30  CONTINUE 

Check  if  In-Plane  Thrust  Exceeds  Drag  (over  complete  orbit)  for 
the  Input  Parameters  and  Raise  Initial  Orbit  Altitude  as  Needed 

ALT  =  ALTI 

IX)  WHILE(TEFFTE.0.D0) 

R  =  KTM’ALT  +  RE 
VEL2  =  MU/R 

CALL  DENS76(ALTT)ENSITY) 

DENSITY  =  FDENS*DENSITY/KTM3 

DRAG  »  0.5DO*CD*AREA*DENSITY*VEL2 

Estimate  Fraction  Orbit  is  in  Sunlight  for  Worst-Case  Date 

FS  *  0.5D0  +  DACOS(RE/R)/PI 

TEFF  =  CALPH*THRUSTI  -  DRAG/FS 

Compute  Net  Effective  Thrust  with  10%  Minimum  Margin  over  Drag 
TEFF  *  0.9D0*CALPH*THRUSTI  -  DRAG/FS 
F(TEFFXT.0.D0)  ALT  =  ALT  +  10.D0 
END  DO 

IF(ALT.GT.ALTI)  THEN 
WRTTE(6,40)  ALTI  ALT 

4  0  FORMAT(/,48H*  DRAG  EXCEEDED  90%  OF  IN-PLANE  THRUST  FOR  INPUT, 

*  42H  PARAMETERS;  STARTING  ALTITUDE  RAISED  FROM.F9.3.6H  KM  TO, 

*  F9.3.3H  KMJ) 

ALTI  =  ALT 

END  IF 

Perform  Spiral  Orbit  Transfer  Analysis 

CALL  TRANSFER 
CLOSE  (UNIT=6) 

IF(IPLOTEQ.  1)  CLOSE  (UNIT=8) 

END 
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SUBROUTINE  TRANSFER 
C 

C  Subroutine  TRANSFER  determines  time  and  delta-V  requirements  for  a 
C  low-thrust  spiral  transfer  from  low  earth  orbit.  Includes  effects 
C  of  atmospheric  drag,  earth  shadow  eclipses,  solar  array  degradation, 

C  and  inclination  change.  Tangential  in-plane  thrusting  and  constant 
C  out-of-plane  thrust  angle  are  assumed  during  transfer.  The  transfer 
C  orbit  is  also  assumed  to  remain  nearly  circular  at  any  instant. 

C 

IMPLICIT  DOUBLE  PRECISION  (A-M.O-Z) 

INTEGER  IDRAGJDEGRD,ICELL,IST(2),IPLOT,IPRINT,IFL,IAL,ICL 
INTEGER  ITER.ILAST 
DIMENSION  FLUEN(5),FST(2) 

COMMON/CONST/MU, SQMUJIE,PIPI2,DTS,DTR,KTM,KTM3,G0,DYEAR,SEPSLN 
COMMON/FLAGS/IDRAG,IDEGRD,ICELL,IST,FST,IPRINT,IPLOT,NPRINT 
COMMON/TRAJ/T.DAYS, DATE  ,ALT,ALTI,ALTF,VEL,R  DOT, INCL.INCLIJNCLF, 

*  ALTINCL,CALPH,OMEGAI,BETA,FE,CD,AREA,FDENS,GAMMA,ORBITS,PERIOD, 

*  FPWR,FLUEN,CFLUEN,THRUST,THRUSTI,THRUSTP,DRAG,ISP,ETA,DVTOT, 

*  DRGLOSS,GRVLOSS,MASSI,MASS,MPROP,SYSPWR,EPPWR,POWER,TECLPS, 

*  NUMTHST 

COMMON/MASS/STR,THMGMT,INEFF,PROPSYSJEPPWER,SASM,MEOTVMPAY,AAREA, 

*  AMASS,XSPWR,PWRWGT,PWRDEN,SCMASS,PPUEFF,ENGEFF,STRFRACT 
COMMON/S  UN/TS  UN 
DATA  ILAST,RE2/1,4.0680606D14/ 

Set  Parameters  for  Inclination  Iteration 

ITER  =  0 

SALPH  =  DSQRT(1.D0  -  CALPH*CALPH) 

Set  Flag  to  Inhibit  Output  Until  Final  Iteration  Step 
IF(INCLF.NE.INCLI)  ILAST  =  0 
Set  Initial  Point  Values  for  Secant  Iteration  Scheme 
SALPHO  =  0.D0 
IDELTO  =  INCLI  -  INCLF 


Initialize  Trajectory  Variables  at  Start  of  Transfer 

10  ALT  =  ALTI 

R  =  KTM*ALT  +  RE 
T  =  0.D0 
DAYS  =  0.D0 
FPWR  =  1.D0 
INCL  =  INCLI 
VEL2  =  MU/R 
VEL  =  DSQRT(VEL2) 

MASS  =  MASSI 
MPROP  =  0.D0 
DVTOT  =  0.D0 
OMEGA  =  OMEGAI 
POWER  =  SYSPWR 
THRUST  =  THRUSTI 
ORBITS  =  0.D0 
CFLUEN  =  0.D0 
TSUN  =  0.0D0 
TECLPS  =  0.D0 
DRGLOSS  =  0.D0 
GRVLOSS  =  0.D0 
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C  Set  Interpolation  and  Print  Indices 
IFL  =  2 
IAL  =  2 
ICL  =  3 
NPT  =  1 
NPR  =  0 
NFL  =  0 

Compute  Drag  Force  Acting  on  Spacecraft 

2  0  IF((ALT.LT.l.D3).AND.(IDRAG.EQ.l))  THEN 
ALT  =  (R  -  RE)/KTM 
CALL  DENS76(ALT,DENSITY) 

DENSITY  =  FDENS*DENSITY/KTM3 
DRAG  =  CD*AREA*0.5D0*DENSITY*VEL2 
ELSE 

DRAG  =  0.0D0 
END  IF 

Compute  Fraction  of  Orbit  that  Spacecraft  is  in  Earth’s  Shadow 

Compute  Ascending  Node  Regression  Rate  (radians  per  day) 

O  MEG  DOT  =  -0. 173903DO*(RE/R)**3.5DO*DCOS(INCL) 

Compute  Solar  Declination  and  Right  Ascension  Angles  for  Current  Date 
Cl  =  PI2*(DATE  +  DAYS  -  80.D0)/DYEAR 
DELTA  =  DASIN(SEPSLN*  DSIN(C  1 )) 

RATS  =  D  ACOS(DCOS(C  1  )/DCOS(DELT  A)) 

IF(Cl.GT.PI)  RATS  =  -l.D0*RATS 
Compute  Beta  Angle  Between  Orbital  Plane  and  Sun  Vector 
BETA  =  DASIN(DCOS(DELTA)*DSIN(INCL)*DSIN(OMEGA  -  RATS) 

*  +  DSIN(DELTA)*DCOS(INCL)) 

IF(DABS(BETA).LT.DASIN(RE/R))  THEN 
ALTM  =  ALT*1000.D0 

FE  =  DACOS(DSQRT(ALTM*ALTM  +  2.DO*RE*ALTM)/(R*DCOS(BETA)))/P1 
ELSE 

FE  =  0.0D0 
END  IF 

Compute  In-Plane  Thrust  Component  while  Changing  Inclination 

IF((INCLF.NE.INCLI).AND.(ALT.GE.ALTINCL))  THEN 
THRUSTP  =  CALPH*THRUST 
THRUSTN  =  SALPH*THRUST 
ELSE 

THRUSTP  =  THRUST 
THRUSTN  =  0.D0 
END  IF 

Compute  Average  Radius  Rate  of  Change 

R3MU  =  2.0D0*DSQRT(R*R*R/MU) 

RDOT  =  R3MU*(THRUSTP*(1.0D0  -  FE)  -  DRAG)/MASS 

Compute  Orbital  Period  with  Spiral  Orbit  Adjustment  if  RDOT  >  0.1  m/s 

IF(DABS(RDOT).GT.0. 1  DO)  THEN 

PERIOD  =  (R/RDOT)*(MU/(SQMU  -  PI*RDOT*DSQRT(R))**2.DO  -  EDO) 
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ELSE 

PERIOD  =  2.0D0*PI*R**1.5D0/SQMU 
END  IF 

Compute  Altitude  Increase  over  Orbit  and  Check  Against  Final  Altitude 

DALT  =  RDOT*PERIOD/KTM 
ALT  =  ALT  +  DALT 
IF(ALT.GT.ALTF)  THEN 

Adjust  Period  and  Altitude  to  Match  Desired  Final  Altitude 
Cl  =  KTM*(ALT  -  ALTF)/RDOT 
ORBITS  =  ORBITS  +  1.D0  -  C  1/PERIOD 
PERIOD  =  PERIOD  -  Cl 
DALT  =  DALT  +  ALTF  -  ALT 
ALT  =  ALTF 
ELSE 

ORBITS  =  ORBITS  +  1.D0 
END  IF 

Update  Orbital  Parameters  and  Compute  Total  delta-V 

T  =  T  + PERIOD 

DAYS  =  DAYS  +  PERIOD/DTS 

R  =  KTM*ALT  +  RE 

VEL20LD  =  VEL2 

VEL2  =  MU/R 

VELOLD  =  VEL 

VEL  =  DSQRT(VEL2) 

Compute  Mean  Flight  Path  Angle  with  respect  to  Local  Horizontal 
GAMMA  =  2.DO*(R*R/MU)*(THRUSTP*(1.0DO  -  FE)  -  DRAG)/MASS 
Compute  Propellant  Mass  Expelled  during  Previous  Orbit 
DMASS  =  PERIOD*(  1.0D0  -  FE)*THRUST/(GO*ISP) 

Compute  delta-V  for  Previous  Orbit  and  Cumulative  delta-V 
DV  =  GO*ISP*DLOG(MASS/(MASS  -  DMASS)) 

DVTOT  =  DVTOT  +  DV 

Compute  Cumulative  Eclipse  Time  during  Transfer 

TECLPS  =  TECLPS  +  FE*PER10D 

TSUN  =  TSUN  +  (l.DO-FE)*PERIOD 

Compute  Total  Propellant  Expended  and  Current  Spacecraft  Mass 

MPROP  =  MPROP  +  DMASS 

MASS  =  MASS  -  DMASS 

Compute  Shift  in  Orbit's  Ascending  Node  (radians) 

OMEGA  =  OMEGA  +  OMEGDOT*PERIOD/DTS 

Estimate  Cumulative  delta-V  Drag  and  Gravity  Losses 

Compute  delta-V  Needed  for  Same  Orbit  Change  with  No  Drag 
IF(DR AG.GT.O.  DO)  THEN 
THRUSTD  =  THRUSTP  -  DRAG/(1.0D0  -  FE) 

THRUSTO  =  DSQRT(THRUSTN*THRUSTN  +  THRUSTD*THRUSTD) 
DMAS  SO  =  PERIOD*(  1.0D0  -  FE)*THRUSTO/(GO*ISP) 

DVO  =  GO*ISP*DLOG(MASS/(MASS  -  DMASSO)) 

ELSE 

DVO  =  DV 
END  IF 

DRGLOSS  =  DRGLOSS  +  DV  -  DVO 

GRVLOSS  =  GRVLOSS  +  GO*RE2/(R*R)*GAMMA*PERIOD 


29 


nnn  n  n  non  on  non  non  non 


Compute  Inclination  Change  based  on  Relation  in  Edelbaum’s  Paper 

IF((INCLF.NE.INCLI).AND.(ALT.GE.ALTINCL))  THEN 
Cl  =  2.  DO/PI 

DINCL  =  C1*DAC0S((VEL20LD  +  VEL2  -  DVO*DVO)/(2.DO*VEL*VELOLD)) 
INCL  =  INCL  -  DINCL 
ELSE 

INCL  =  INCL 
END  IF 

Compute  Cumulative  Electron/Proton  Fluence 

IF(IDEGRD.GE.  1 )  THEN 
ALT1  =  ALT  -  0.5D0*DALT 

CALL  FLUENCE(ALT1  .INCL, FLUEN.IAL.ICL, 1ST, FST) 

CFLUEN  =  CFLUEN  +  FLUEN(5)*PERIOD/(DTS*DYEAR) 

END  IF 

Determine  Array  Power  and  Thrust  Degradation  due  to  Cum  Fluence 

IF(IDEGRD.EQ.2)  THEN 
CALL  DEG RAD( CFLUEN, FPWR.ICELL.IFL) 

POWER  =  FP  WR  *S  Y  SPWR 
F(POWER.GT.O.DO)  THEN 

(Future  upgrade  to  include  effect  of  Isp  drop  and  logic  for 
possible  step  thrust  decreases  for  multiple  ion  thrusters.) 

THRUST  =  FPWR*THRUSTI 
ELSE 

WRITE(6,30)  ALT 

3  0  FORMAT(//,40H*  ZERO  POWER  AVAILABLE  FOR  THRUSTERS  AT  ,F9.3, 

*  44H  KM  DUE  TO  ARRAY  DEGRADATION;  RUN  TERMINATED,//) 
RETURN 
END  IF 
END  IF 

Output  Data  at  Specified  Time  Steps  during  Transfer 

IF(ILAST.EQ.O)  GOTO  50 
IF(IPRINT.LE.O)  GOTO  40 
NPR  =  NPR  +  1 
IF(NPR.LT.NPRINT)  GOTO  40 
NPR  =  0 

CALL  OUTPUTl(NPT) 

4  0  IF(IPLOT.LE.O)  GOTO  50 

NPL  =  NPL  +  1 
IF(NPL.LT.NPLOT)  GOTO  50 
NPL  =  0 

CALL  PLOT(DAYS,NORBIT,ALT,INCL,DVTOT, CFLUEN, FPWR) 

50  IF(ALT.LT.ALTF)  GOTO  20 
IF(ILAST.EQ.l)  GOTO  70 

Check  Final  Inclination  Convergence  and  Adjust  sin( ALPHA)  Value 
ITER  =  ITER  +  1 
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IDELT  =  INCL  -  INCLF 

.  C  Set  Flag  for  Last  Iteration  Step  if  Inclin  is  within  0.5  deg  of  Goal 
IF(DABS(IDELT).LE.0.0087266463DO)  ILAST  =  1 
C  Estimate  Adjustment  to  sin(ALPHA)  using  Secant  Method 
DALPH  =  IDELT*(SALPH  -  SALPH0)/(IDELT  -  IDELTO) 

SALPHO  =  SALPH 

IDELTO  =  IDELT 

SALPH  =  SALPH  -  DALPH 

CALPH  =  DSQRT(1.D0  -  SALPH*SALPH) 

IF(ITER.LT.10)  GOTO  10 
WRJTE(6,60) 

60  FORMAT(/,5 1 H* SECANT  ITERATION  FOR  FINAL  INCL  FAILED  TO  CONVERGE) 
7  0  IF(IPRINT.EQ.  1 )  CALL  OUTPUT  1  (NPT) 

IF(IPLOT.EQ.l)  CALL  PLOT(DAYS,NORBIT,ALT,INCL,DVTOT,CFLUEN,FPWR) 

C 

AMASS  =  SYSPWR*1000.D0/PWRWGT 
STR  =  STRFRACT  *  (SCMASS+AMASS-MPROP) 

INEFF  =  (1.0D0-PPUEFF)  *  SYSPWR 
SMMEOTV  =  SASM  +  PROPSYS 

MEOTV  =  ((SMMEOTV*EPPWR)+(THMGMT*INEFF))+STR 
MPAY  =  SCMASS+AMASS-MPROP-MEOTV 

Write  Mission  Data  to  Output  File 

WRITE(6f350) 

WRITE(6,400)  THRUS  i  /  ^rLOAT(NUMTHST),NUMTHST,THRUST,ISP 
WRITE(6,450)  PWPP'  N  .PWRWGT, SYSPWR, XSPWR, A  AREA, AMASS.ETA 
IF(IDECAY.GT.O)  Yv  KITE(6,500)  ALTMIN 

WRITE(6,575)  JTR,THMGMT*INEFF,PROPSYS*EPPWR,SASM*EPPWR,MEOTV, 

+  MPA  Y.MEOTV +MPA  Y  *MPROP,SCMASS+ AMASS 

W'RITE(6,550)ALTI,ALTF,DVTOT,TSUN/DTS,TECLPS/DTS, 

+  DAYS,  ORBITS 

format  statements... 

350  FORMAT(//,20x, 'ELECTRIC  PROPULSION  MISSION  SUMMARY') 

400  FORMAT(//,’  THRUSTERS', 

+  //,’  Thrust/Thruster  (N):  ',2x,fl7.2, 

+  /,'  Quantity:  ',5x,i8, 

+  /,' Total  Thrust  (N):  ',2x,fl7.2, 

+  /,’  Specific  Impulse  (sec):',4X,F16.0) 

450  FORMAT(/,' POWER:’, 

+  //,'  Power  Density  (W/M**2):’,2x,f7.2, 

+  /,'  Specific  Power  (W/kg):  ’,2x,f7.2, 

+  I,'  Total  Power  (kW):  ’,2x,n.2, 

+  /,’  Excess  Power  (kW):  ',4x,f5.2, 

+  /,’  Array  Area  (M**2):  ’,lx,f8.2, 

•f  /,' Array  Mass  (kg):  ',lx,f8.2, 

+  /,’  System  Efficiency:  \5X,F5.3) 

500  FORMATS, ’  ORBIT  TRANSFER  WILL  BEGIN  AT\lx,fl2.2,'(km)',/, 

+  IN  ORDER  TO  OVERCOME  THE  DRAG  FORCE  ) 

550  FORMAT(/,’ MISSION:', 

+  //,’  Initial  Altitude  (km):  ',3x,f7.0, 

+  /,’  Final  altitude  (km):  ’,3x,f7.0, 

+  /,'  Required  Delta-V  (m/s):',2x,f7.2, 

+  //,’  Time  in  sun  (days):  ',2x,f7.2, 

+  /,'  Time  in  shadow  (days):  ',2x,f7.2, 
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+  /,'  Total  trip  time(days):  ’,2x,F7.2 

+  /,’  Total  number  of  orbits:’, 2x, FI 0.2) 

575  FORMAT(/,'  SPACECRAFT  WEIGHT  BREAKDOWN  (kg):’, 
+  //,'  Structural  mass:  ’,lx,f7.2, 

+  /,'  Thermal  management  mass:',lx,f7.2, 

+  /,’  Propulsion  system  mass:  ’,lx,f7.2, 

+  /,’ Power  mass:  \lx,f7.2, 

+  /,19x,f8.2,’  *Total  EOTV  dry  mass’, 

+  /,’  Payload:  \lx,f8.2, 

+  /,19x,f8.2,’  **Total  Mass  @  GEO’, 

+  /,’  Propellant:  \lx,f8.2, 

+  /,19x,f8.2,’  ***Total  Mass  @  LEO’) 

600  FORMAT(/,’  Time  penalty  due  to  drag:’,lx,f7.2,’  days’, 

+  /,’  Delta  V  penalty  due  to  drag:’,lX,F7.2,’  m/sec’) 

C 

RETURN 

END 
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SUBROUTINE  DEGRAD(CFLUEN, FPWR, ICELL, IFL) 

.  C 

C  Subroutine  DEGRAD  computes  solar  array  power  degradation  based  on  the  cumulative  electron  and 
C  proton  fluence  experienced  during  the  orbit  transfer.  Normalized  cell  power  data  stored  for 

C  cumulative  1  MeV  electron  fluences  of  1.E12,  3.E12,  1.E13,  3.E13,  ...  ,  1.E16.  Fluence  values 

C  stored  as  logarithms  to  facilitate  interpolatioa 
C 

IMPLICIT  DOUBLE  PRECISION  (A-M.O-Z) 

INTEGER  I.IFL, ICELL 

DIMENSION  LFLUEN(9),  PWRNORM(9,8) 

C  Reference  Cumulative  Fluences  (Natural  Logarithm  Values)  for  which 
C  Solar  Cell  Normalized  Power  Data  is  Input  -  IE  12  to  IE  16  MeV 

DATA  (LFLUEN(I),  I  =  1,9)  /27.631021D0,  28.729633D0,  29.933606D0, 

*  31.032218DO,  32.236191DO,  33.348043D0,  34.538775DO,  35.637389D0, 

*  36. 84 136 IDO/ 

C  Normalized  Max  Power  vs.  Fluence  for  Si  BBSFR  Thin  Cell,  3.4  mils; 

C  "Solar  Cell  Radiation  Handbook",  Addendum  1,  15  Feb  89,  Figure  18. 

DATA  (PWRNORM(I,l),  I  =  1,9)  /1.000D0,  0.992D0,  0.97 IDO,  0.938D0, 

*  0.889D0,  0.821D0,  0.726D0,  0.632D0,  O.530DO/ 

C  Normalized  Max  Power  vs.  Fluence  for  ASEC  OMCVD  GaAs/Ge  Cell; 

C  "Characterization  of  GaAs  Solar  Cells",  JPL  Pub.  88-39,  Fig.  10. 

DATA  (PWRNORM(I,2),  I  =  1,9)  /1.000D0,  0.994D0,  0.976D0,  0.957D0, 

*  0.923D0,  0.874D0,  0.801D0,  0.691D0,  0.521D0/ 

IF(CFLUEN.LE.  1  .D 1 2)  GOTO  30 


Interpolate  for  Normalized  Power  Degradation  Factor 


LCFLUEN  =  DLOG(CFLUEN) 

1  0  IF((LCFLUEN.LE.LFLUEN(IFL)).OR.(IFL.EQ.8))  GOTO  20 
IFL  =  IFL  +  1 
GOTO  10 

20  Cl  =  (LCFLUEN  -  LFLUEN(IFL- 1  ))/(LFLUEN(IFL)  -  LFLUEN(IFL-l)) 

FPWR  =  PWRNORM(IFL-l.ICELL)  +  Cl*(PWRNORM(IFL, ICELL) 

*  -  PWRNORM(IFL- 1  ,ICELL)) 

RETURN 

30  FPWR  =  1.0D0 
RETURN 
END 
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SUBROUTINE  FLUENCE(ALT,INCL,FLUEN,IAL,ICL,IST,FST) 

C 

C  Subroutine  FLUENCE  computes  annual  equivalent  1  MeV  Electron  Fluence  from  Trapped 

C  Electrons  (Pm ax  and  Voc)  and  Trapped  Protons  (Pm ax  and  Voc)  as  a  function  of  spacecraft  orbital  C 

altitude,  inclination,  and  equivalent  array  front  and  back-side  shield  thicknesses.  (Note:  Fluence 
C  data  covers  altitudes  of  277  to  35,794  km,  inclinations  of  0  to  30  degrees,  and  0  to  60  mils  shield 

C  thickness.  Fluence  data  for  higher  inclination  orbits  maybe  added  by  raising  inclination  index 

C  beyond  4  and  expanding  tables.) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-M.O-Z) 

INTEGER  I,I1,I2,IS,IAL,ICL,IST(2) 

DIMENSION  ALTR(34),INCLR(4),FLUEN(5),FST(2) 

DIMENSION  EMAX(34,4,8),PMAX(34,4,8) 

C  Reference  Inclinations  for  Fluence  Data,  radians 
DATA  (INCLR(I),  I  =  1,4) 

*  /  0.0D0,  0.1745329252D0,  0.3490658504D0,  0.5235987756D0/ 

C  Reference  Altitudes  for  Fluence  Data,  km 

DATA  (ALTR(I),  I  =  1,34) 

*  /  277D0,  463D0,  555D0,  833DO,  1111D0,  1481D0,  1852D0, 

*  2315D0,  2778D0,  3241D0.  3704D0,  4167D0,  4630D0,  5093D0, 

*  5556D0,  6482D0,  7408D0,  8334D0,  9260D0,  101 85D0,  1 1 1 12D0, 

*  12964 EX),  14816DO,  16668D0,  18520DO,  20372D0,  22224D0,  24076D0, 

*  25928D0,  27780D0,  29632D0,  31484D0,  33336D0,  35794D0/ 

C  Annual  Equivalent  1  MEV  Electron  Fluence  from  Trapped  Electrons  for 
C  0  mils  Shield  Thickness  and  Orbit  Inclinations  of  0,  10,  20,  30  deg 
DATA  ((EMAX(I1,I2,1),  II  =  1,34),  12  =  1,4) 

*  /0.00D00,  1.14D07,  4.65D07,  8.03D09,  4.92D1 1,  7.25D12,  2.49D13, 

*  6.03D13,  9.80D13,  1.26D14,  1.46D14,  1.64D14,  1.73D14,  1.70D14,  ... 


(3  pages  of  numerical  tables  have  been  omitted  here.  See  reference  for  a  complete  listing  of  these  tables) 
Find  Interpolation  Coefficients  for  Increasing  Alt/Decreasing  Inclin 

1  0  IF((ALT.LE.  ALTR(IAL)).OR.(I AL.  EQ.  34))  GOTO  20 
IAL  =  IAL  +  1 
GOTO  10 

2  0  FALT  =  (ALT  -  ALTR(IAL- 1  ))/(ALTR(I AL)  -  ALTR(IAL- 1 )) 
IF((INCL.GE.INCLR(ICL)).OR.(ICL.EQ.  1))  GOTO  30 
ICL  =  ICL  -  1 
GOTO  20 

3  0  FIN CL  =  (INiJL  -  1NCLR(ICL))/(INCLR(ICL+1)  -  INCLR(ICL)) 

Loop  to  Compute  Fluences  on  Front  (IS=1)  and  Back-Side  (IS=2)  of  Array 

DO  50  IS  =  1,2 

C  Compute  Trapped  Electron  Fluence  (Pmax,  Voc,  Isc)  as  Function  of 
C  Orbital  Inclination,  Altitude,  and  Solar  Array  Shield  Thickness 
I  =  IST(IS) 
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CAl  =  EMAX(IAL- 1  ,ICL,I)  +  FINCL*(EMAX(IAL-1,ICL+1,I) 

*  -EMAX(IAL-UCL.I)) 
CA2  =  EMAX(IAL,ICL,I)  +  FINCL*(EMAX(IAL,ICL+1,I) 

*  -  EMAX(IAL,ICL,I)) 
FLUEN(IS)  =  CAl  +  FALT*(CA2  -  CAl) 

C  Bypass  Interpolation  Calculations  if  Ref.  Shield  Thickness  Value  Used 
IF(ABS(FST(IS)).LT.0.001)  GOTO  40 
I  =  IST(IS)  +  1 

CAl  =  EMAX(I  AL- 1  ,ICL,I)  +  FINCL*(EMAX(IAL-1,ICL+1,I) 

*  -  EMAX(IAL-1,ICL,I)) 
CA2  =  EMAX(IAL,ICL,I)  +  FINCL*(EMAX(IAL-1,ICL+1,I) 

*  -  EMAX(IAL.ICL.I)) 
FLUEN(IS)  =  FLUEN(IS)  +  FST(IS)*((CA1  +  FALT*(CA2  -  CAl)) 

*  -  FLUEN(IS)) 


Compute  Trapped  Proton  (Pmax  and  Voc)  Fluence  as  a  Function  of 
Orbital  Inclination,  Altitude,  and  Solar  Array  Shield  Thickness 

40  I  =  IST(IS) 

12  =  IS  +  2 

CAl  =  PMAX(I  AL- 1  ,ICL,I)  +  FINCL*(PMAX(IAL-1,ICL+1,I) 

*  -  PMAX(IAL-1,ICL,I)) 
CA2  =  PMAX(IAL.ICL.I)  +  FINCL*(PMAX(I  AL.ICL+ 1 ,1) 

*  -  PMAX(IAL,ICL,I)) 
FLUEN(I2)  =  CAl  +  FALT*(CA2  -  CAl) 

Bypass  Interpolation  Calculations  if  Ref.  Shield  Thickness  Value  Used 
IF(ABS(FST(IS)).LT.0.001)  GOTO  50 
I  =  IST(IS)  +  1 

CAl  =  EMAX(I  AL- 1  ,ICL,I)  +  FINCL*(EMAX(IAL-1,ICL+1,I) 

*  -  EMAX(IAL-1,ICL,I)) 
CA2  =  EMAX(I  AL.ICL.I)  +  FINCL*(EMAX(I  AL.ICL+ 1 ,1) 

*  -  EMAX(I AL.ICL.I)) 
FLUEN(I2)  =  FLUEN(I2)  +  FST(IS)*((CA1  +  FALT*(CA2  -  CAl)) 

*  -  FLUEN(I2)) 
50  CONTINUE 

Compute  Combined  Fluence  on  Front  and  Back  Sides  of  Array 

FLUEN(5)  =  FLUEN(l)  +  FLUEN(2)  +  FLUEN(3)  +  FLUEN(4) 

RETURN 

END 
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SUBROUTINE  PLOT(DAYS,NORBIT,ALT,INCL,DVTOT,CFLUEN,FPWR) 
IMPLICIT  DOUBLE  PRECISION  (A-M.O-Z) 

DATA  DTR/0.17453292519943D0/ 


Write  Trajectory  Data  to  Plot  File  8 
INCLD  =  INCL/DTR 

WRTTE(8, 100)  DAYS, ORB  ITS, ALT, INCLD, DVTOT  FPWR.CFLUEN 
100  FORMAT(6F12.4,D12.4) 

RETURN 

END 
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